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, Abstract 

In the resonance model, high-frequency quasi-periodic oscillations (QPOs) are supposed to be 
a consequence of nonlinear resonance between modes of oscillations occurring within the innermost 
parts of an accretion disk. Several models with a prescribed mode-mode interaction were proposed 
in order to explain the characteristic properties of the resonance in QPO sources. In this paper, we 
examine nonlinear oscillations of a system having a quadratic nonlinearity and we show that this 
, case is particularly relevant for QPOs. We present a very convenient way how to study internal 

00 ■ resonances of a fully general system using the method of multiple scales. Finally, we concentrate 

' to conservative systems and discuss their behavior near the 3 : 2 resonance. 

o 

1 Introduction 

Of 

there is a natural and attractive possibility of explaining the observed rational ratios of high-frequency 
QPOs as a consequence of non-linear coupling between di fferent modes of accretion disk oscillatio ns. 
The idea has been pursued in several papers (recently, e.g. lAbramowicz et al.ll20o5 lRebuscoll2004j) . 

Specific models invoke particular physical mechanisms. Some models can be almost immediately 
comprehended as distinct realisations o f the general approach discussed here - for example, various 
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In the resonance model l|Abramowicz fc Kliiznialdl200ll iKlnzniak fc Ahra,mowic7ll2000l: iKatol (20031. 



formulations of the orbiting spot model IjSchnittman fc Bertchingerll2004l) or the models, where QP Os 



are produced by the magnetically driven resonance in a diamagnetic accretion disk l|Lailll999|) 
while other seem to be more distant from the view presented herein - e.g. the transiti on layer model 
llTitarchukll2002|) . an interesting idea of p-mode oscill ations of a small accretion torus ijRezzola et alJ 
2003) or the model of blobs in an accr etion d i sc (se e e.g. lKarasll999l:lLi fc Naravanl2004 and references 
cited therein). Also in this context. IKatol l|2004f) discussed the resonant interaction between waves 
propagating in a warped disk, including their rigorous mathematical description. Instead of pursuing 
a specific model, here we keep the discussion as general as possible, aiming to implement the formalism 
of multiple scales. Indeed, we show that there is unquestionable appeal in this approach which offers 
some additional insight into generic properties of resonant oscillations. 

Some properties of an accretion disk oscillations can be discussed within the epicyclic approxima- 
tion of a test particle on a circular orbit near equatorial plane. Suppose that angular momentum of 
the particle is fixed to a value £. The effective potential Ue(r,9) has a minimum at radius rrj, corre- 
sponding to the location of the stable circular orbit. An observer moving along this orbit measures 
radial, vertical and azimuthal epicyclic oscillations of a particle nearby. Since the angular momentum 
of the particle is conserved, only two of them - radial and vertical - are independent. The epicyclic 
frequencies can be derived from the geodesic equations expanded to the linear order in deviations 

*To appear in proceedings of the Workshop Processes in the vicinity of black holes and neutron stars, held at Silesian 
University (Opava, October 2003). 
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5r = r — r*o and 89 = 9 — it/ 2 from the circular orbit. We get two independent second-order differential 
equations describing two uncoupled oscillators with frequencies u> r and ujg, which are given by the 
second derivatives of effective potential Ut(r, 9). In Newtonian theory, u> r and ojg are equal to the Ke- 
plerian orbital frequency £Ik- This is in tune with the fact that orbits of particles are planar and closed 
curves. The degeneracy between two ep icyclic frequencies can be seen as a result of scale-freedom 
of the Newtonian gravitational potential (Ab ramowicz & Kluzn iak 2003). In Schwarzschild geometry 
this freedom is broken by introducing the gravitational radius r g = 2GM/c 2 . The degeneracy between 
the vertical epicyclic and the orbital frequencies is related to spherical symmetry of the gravitational 
potential, which assures the existence of planar trajectories of particles. All three frequencies are 
different in the vicinity of a rotating Kerr black hole. 

In addition, when nonlinear terms of geodesic equations are included, the two oscillations in r 
and 9 directions become coupled and variety of new phenomena connected to nonlinear nature of 
the equations appear. This rich phenomenology includes frequency shift of observed frequencies with 
respect to eigenfrequencies, presence of higher harmonics and subharmonics, drifts and parametric 
resonance. The first three are connected to nonlinear oscillations of each mode and the last one comes 
from the coupling between two modes. 

The paper is organized in a following way: In the next section we will study one-deg ree-of- freedom 
system with a quadra tic nonlinearity. We introduce the method of multiple scales l|Navfehlll973t 
iNavfeh fc Mooklll979l) and, using this formalism, we derive key properties of nonlinear oscillations. 
In the third section we perform the multiple-scales expansion in the case of a fully general system 
of two coupled oscillators and find possible parametric resonances up to the fourth order. Then we 
concentrate on the 3 : 2 parametric resonance of the conservative system. 



2 Effects of nonlinearities 

Let us consider the case of small but finite oscillations of a single-degree-of-freedom system with 
quadratic nonlinearity governed by equation 

•■ i 2 2 2 /-,\ 

x + oj x = aoj x . (1) 

The strength of the nonlinearity is parametrized by constant a. When a = one obtains governing 
equation of the corresponding linear system. We seek a perturbation expansion of the form 

x(t, e) = exi(t) + e 2 x 2 (t) + e 3 x 3 (t) + C(e 4 ), (2) 

The expansion parameter e expresses the order of amplitude of oscillations. The main advantage of 
this approach is that, although the original equation is nonlinear, we solve linear equations in each 
step. For a practical purpose we require this expansion to be uniformly convergent for all times of 
interest. In that case the higher-order terms are small compared to lower-order terms and a sufficient 
approximation is reached concerning the finite number of terms. The expansion J2J can represent 
a periodic solution as well as an unbounded solution with exponential grow. The uniformity of the 
expansion means that the higher-order terms are not larger than the corresponding lower-order ones. 

We substitute expansion (J2J) into governing equation JQ) and, since Xk is independent of e, we 
equate coefficients of corresponding powers of e on both sides. This leads to the following system of 
equations 

Xl + lo 2 xi = 0, (3) 
X2+lo 2 x 2 = au 2 x\, (4) 
X3+oj 2 X3 — 2a<jj 2 xiX2- (5) 

The general solution of eq. © can be written in the form xi(t) = Ae luJt +cc, where cc denotes complex 
conjugation. The complex constant A contains information about the initial amplitude and phase of 
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oscillations. Substituting it into eq. Q we find linear equation for the first approximation x 2 (i) 

x 2 +u; 2 x 2 =a(A 2 e 2 ^ t + \A\ 2 )+cc. (6) 
A general solution consists of the solution of the homogeneous equation and a particular solution, 

x 2 (t) = A 2 e^ - a Q^e*** - \A X \ 2 \ + cc, (7) 

where A\ denotes a constant A of the solution of eq. J3J). Therefore, the solution of governing equation 
up to the second order is given by 

x{t) = (eA 1 + e 2 A 2 ) - a Q^e 3 *"* - |^i | 2 ) + cc. (8) 

In fact, there are two possible ways how to satisfy general initial conditions a;(0) = exo and i(0) = eio 
imposed on equation Q). The first one is to compare them with the general solution JSJ and find 
constants A±, A 2 . This procedure should be repeated in each order of approximation which involves 
quite complicated algebra especially in higher orders. The second, equivalent and apparently much 
easier way is to include only particular solutions to the higher approximations and treat the constant 
A as a function of e with expansion A = A\ + tA 2 + . . .. Then, given initial conditions are satisfied 
by expanding the solution for X\ via e and choosing the coefficients A n appropriately. 
According to this discussion we express the solution of eq. © as 

x 2 (*) = -a (\A 2 e 21 ^ - | A\ 2 ) + cc. (9) 
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Substituting x\ and x 2 into (JSJ we obtain 



£a + J* X3 = ^! ( 5 A| A\ 2 e wt - A 3 e 31 ^) + cc. (10) 

o 

Since the right-hand side of this equation contains the term proportional to e llAjt , any solution must 
contain a secular term proportional to te lut , which becomes unbounded as t — > oo. This fact has 
nothing to do with true physical behavior of the system for large times. The meaning is rather 
mathematical. Starting from time when (tot) ~ l/(ea), the higher-order approximation x$, which 
contains the secular term, does not provide only a small correction to x\ and x 2 , and the expansion 
becomes singular. However, the presence of the secular term in the third order reflects very general 
feature of nonlinear oscillations - dependence of the observed frequency on the actual amplitude. For 
larger amplitudes the actual frequency of oscillations differs from the eigenfrequency to and the higher- 
order terms in the expansion J5J - always oscillating with an integer multiples of w — must quickly 
increase as time grows. 

2.1 The method of multiple scales 

Is it possible to find an expansion representing a solution of equation Q which is uniformly valid even 
for larger time then ~ e _1 ? Yes, if one considers more general form of the expansion than eq. (0). 
In the method of multiple scales more general dependence of coefficients Xi on the time is reached by 
introducing several time scales T^, instead of one physical time t. The time scales are introduced as 

T^ = e% M = 0,l,2... (11) 

and they are treated as independent. It follows that instead of the single time derivative we have an 
expansion of partial derivatives with respect to the 

j t = D a + eD 1 +e 2 D 2 + e :i D 3 + 0(e 4 ), (12) 

rf 2 

— = D 2 +2eD D 1 + e 2 {D 2 +2D D 2 ) + 2e 3 (D D s +D 1 D 2 ) + O(e 4 ), (13) 
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where — d/dT^. 

We assume that the solution can be represented by an expansion having the form 

x(t,e) = ex 1 (T^+e 2 x 2 (T f ,) + e 3 x 3 (T, l ) + 0(e 4 ). (14) 

The number of time scales is always the same as the order at which the expansion is truncated. Here 
we carry out the expansion to the third order and thus first three scales To, T\ and T2 are sufficient. 

Substituting eqs. ((Hfl and (fHj|l into the governing equation and equating the coefficients of e, 
e 2 and e 3 to zero we obtain 

{D 2 +u 2 ) Xl = 0, (15) 
{Dl+u 2 )x 2 = -2D D lXl +auj 2 X!, (16) 
(Dl+uj 2 )x 3 = -2D D x x 2 - D\xx - 2D D 2 x 1 + auj 2 X!X 2 . (17) 

Note that although these equations are more complicated than they are still linear and can be 

solved successively. The solution of the equation lfT5|l is the same as the solution of the corresponding 
linear system, the only difference is that constant A now generally depends on other scales 

xi = A(T 1 ,T 2 )e lujT ° +cc. (18) 

Substituting x% into the equation (|16|) we obtain 

{Dl + lu 2 )x 2 = -2iuj(D 1 A)e iu)T ° + auo 2 [A 2 e 2wTa + \A\ 2 ) + cc. (19) 

The first term on the right-hand side implies the presence of a secular term in the second-order 
approximation which causes non-uniformity of the expansion (|14[) . However, in case of the method 
of multiple scales these terms can be eliminated by imposing additional conditions on the function 
A(T^). These conditions are sometimes called conditions of solvability or conditions of consistency. 
In fact, the reason why the same number of scales as the order of approximation is needed is that one 
secular term gets eliminated in each order, and therefore the function A(T^) is specified by the same 
number of solvability conditions as the number of its variables. 

The secular term is eliminated if we require D\A = 0. In further discussion we assume that A is 
a function of T2 only. A particular solution of equation I jlfljl is 

x 2 (t) = -a [^A 2 (T2)e 2 ^ To - \A{T 2 )\ 2 ^j + cc. (20) 

Using the condition D\A = the equation Ijl7|l takes much simpler form 

(D 2 + u 2 )x z =-( 2zujD 2 A + A\A\ 2 ) e« uT » - ^! A 3 e 3i " T « + cc. (21) 
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The secular term is eliminated equating the terms in the bracket to zero 

2iuD 2 A+ Wa2u}2 A\A\ 2 = 0. (22) 
3 

This additional condition fully determines (excepting initial conditions) time behavior of function 
A(T 2 ). Let us write it in the polar form A = \ae l(f ' and then separate real and imaginary parts. We 
obtain 

5a 2 

D 2 a = and D 2 (p = — -ua 2 . (23) 

1 z 



The solutions of these equations are 
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a and cj) = -— a 2 alujT 2 + <f>o, (24) 
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Figure 1: Oscillations of the system with quadratic nonlinearity governed by the equation 0J. We 
compare results of the multiple scales method (solid curve), a simple straightforward expansion method 
(dashed curve) and the direct numerical integration of the governing equation (points). The initial 
condition is x(0) = 0.1, ±(0) = and the strength of nonlinearity is a = 3. The horizontal dotted line 
shows the shifted 'equilibrium position' and the vertical one denotes the value (cot) = (ea)" 1 at which 
the straightforward expansion becomes nonuniform. The solution corresponding to the linear system 
is also shown (dotted curve). 



where do and </>o are constants which are determined from the initial conditions. 

It follows from eq. 118fl that A(T2) slowly modulates the amplitude and the phase of oscillations. 
Since a is constant, the amplitude is constant all the time. Since depends on Ti = e 2 t linearly, also 
the observed frequency of the oscillations is constant, but not equal to the eigenfrcqucncy to. 

Substituting eqs. and lj2U|l into ijTijl . we obtain solution of eq. (Q) up to the second order 

x(t) = a cos{to*t + (f> ) - %al cos[2(u*i + O )1 + ^Qn + O(an), (25) 

D 2 

where a = ea <C 1 and to* is the observed frequency of oscillation given by 

w* = u;(l-^a§V (26) 

2.2 Essential properties of nonlinear oscillations 

We close this section summarizing main properties of nonlinear oscillations. The equation l|25|) is very 
helpful for this purpose. 

The leading term of the expansion 125JI describes oscillations with frequency close to eigenfrequency 
of the system. Both the amplitude and the frequency are constant in time, but they are not indepen- 
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dent (as in case of the linear approximation). The frequency correction given by (|26[) is proportional 
to the square of the amplitude. This fact is sometimes called the amplitude-frequency interaction and 
- as was mentioned above - it causes the non-uniformity of the expansion @ (see the figure 

The second term oscillates with twice as large frequency and provides the second-order correction 
to the leading term. The presence of higher harmonics is another particular feature of nonlinear 
oscillations and the fact that it has been reported in several sources of QPOs, points to nonlinear 
nature of this phenomena. 

The third term describes constant shift from the equilibrium position and is related to the asym- 
metry of the potential energy about the point x = 0. In the linear analysis, this effect is not present 
because the potential energy depends on x 2 , and therefore is symmetric with respect to the point 
x = 0. Hence, the drift is the third characteristic feature of nonlinear oscillations. 



3 Parametric resonances in conservative systems 

Let us study nonlinear oscillations of the system having two degrees of freedom, i.e., the coordinate 
perturbations 8r and 86. The oscillations are described by two coupled differential equations of the 
very general form 

8r + uj 2 8r = lo 2 f r (8r,89,6r,86), (27) 
56 + u 2 86 = to 2 f e (6r, 66, 8r, 86). (28) 

Suppose that the functions f r and fg are nonlinear, i.e., their Taylor expansions start in the second 
order. Another assumption is that these functions are invariant under reflection of time (i.e., the 
Taylor expansion does not contain odd powers of time derivatives of 8r and 86). As we see later, this 
assumption is related to the conservation of e nergy in the system. M any authors studied such systems 
with a particular form of functions / and g ((Navfeh fc Mook| ^979). however, in this paper we keep 
discussion fully general. 

We seek the solutions of the governing equations in the form of the multiple-scales expansions 

4 4 

8r(t, e) = J2 e n r n (TJ, 86{t, e) = £ t n 6 n {T^), (29) 

n— 1 n— 1 

where = e^t are independent time scales, /i = 0, 1, 2, 3 (we finish the discussion in the fourth order, 
however, it is possible to proceed to higher orders in suggested way). We expand the time derivatives 
according to eqs. (|12|l and (| 1 3|) and equate terms of the same order in e on both sides of the governing 
equations. 

In the first order we obtain equations corresponding to the linear approximation 

(D 2 + uj 2 ) n = 0, (D 2 + oj 2 )6x = 0. (30) 

with the solutions 

xx = A r + A- r , 6x=A e + A- e , (31) 

where we denoted A x = A x e luJxT ° , A- x = A*e~ luJ:cTo , and x = r, 6, respectively (since many consider- 
ations are independent of the mode of oscillations, we keep this notation through the whole section). 
The complex functions A x generally depend on the time-scales T±, T 2 , T 3 . 

Having solved the linear approximation, we can proceed to higher orders. The terms proportional 
to e 2 in the expanded left-hand side of the governing equation H27J) resp. I|28|l are 



Sx + uj 2 x 



2 = i D o+ u 2 x )x 2 + 2iu x DxA x - 2iuj x DxA^ x , (32) 
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tog : u> r Secular terms Solvability condition 











Outside 


D\A r 


D\A r 


= 


resonance 




DxAe 


= 


1 : 2 


A-rAe 


-2iD x A r + u} r KA_ r A e 


= 






-2iD x A e + ujgLAj 


= 


2 : 1 


A 


-2iD x A r + u> r KA$ 


= 




A r A-e 


-2iD x A g + LA r A_ g 


= 



Table 1: Possible resonances and appropriate solvability conditions in the second order of approxima- 
tion. We substitute constants Cq™' 2 ^ by K and L for simplicity. The first record is related to the case 
when the system is far from any listed resonance. In this case only strictly secular terms are present. 
The first resp. second row in the record of each resonance is related to the equation for the radial 
resp. vertical oscillations. In that case, we list only nearly secular terms in the 2nd column, however 
strictly secular terms are included in the solvability conditions. 



On the right-hand side there are second-order terms of the Taylor expansion of the nonlinearity 
f(5r,59,5r,89), with n, 9%, D$r x and D$9 X in the place of Sr, 68, 8r and 69, respectively. Since the 
derivative with respect to To only adds the coefficient ioj x , the second-order terms on the right-hand 
sides can be expressed as linear combinations of quadratic terms constructed from A± r and A±g 

~f x (5r, 69, 8r, 60)] 2 =<4^ A?* A^ A a J r A%, (33) 

M= 2 

where a — (a x , . . . , a^) and \a\ =«!+... + 04. The constants are combinations of u> x and 

the coefficients in the Taylor expansion of functions / and g, respectively. The coefficients coming 
from the terms containing time derivatives are generally complex, since each time derivative produces 
one V . However, if we suppose that the Taylor expansion does not contain odd powers of time 
derivatives, all of the constants Ca' x ^ must be real. Equating right-hand sides of equations (|32|l and 
(13*31 we have 

(Dg + loI)x 2 = -Uu x D x A x + 2ilj x D 1 A_ x +u% £ C^A^A^A1%A%. (34) 

|a|=2 

The right-hand side of equation (|34fl contains one secular term independently of the eigenfrequen- 
cies uj r and ujg. However, in some particular cases, additional secular terms appear. For example, when 
uj r fh 2u>g the terms proportional to Ag in the radial equation (x — > r) and A r A_g in the meridional 
equation (x — > 9) become nearly secular and they should be included in the solvability conditions. 
The analogical situation happens when co r kwj/2. These cases show qualitatively different behavior, 
and we speak about parametric or internal resonance. Possible resonances in the second order of 
approximation and appropriate solvability conditions are listed in the tabled At this moment, let us 
assume that the system is far from the resonance and require 

D X A X = 0. (35) 

In this case the frequencies and the amplitudes are constant and the behavior of the system is almost 
the same as in the linear approximation. The only difference is the presence of the higher harmonics 
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UJg \ UJ r 


Secular terms 


Solvability condition 




Outside 
resonance 


D 2 A r , \Ar\ Ar, \A^\A r 
D 2 A e , \A r \ 2 A e , \A 2 \A e 


2iD 2 Ag -UJg [X r \A r \ 2 - 


i „ U 121 J 
t Kg\Ag\ J A r 

'Xg\Ag\ 2 ' \Ag 


= 
= 


1 : 3 




o;r> 2 1 1 r*- l A I 2 _i_ u-„ I a „ i 2 l 

AlL) 2 A r — UJ r [K r \Ji r \ + Kg\/ig\ J 
2lD 2 Ag - UJg [X r \A r \ 2 + \g\Ag\ 2 ] Ag 


/f j_ / 1 j< 2^3 

/i r -j" UJ r I\ I\q 
4- UJgLA r A 2 _ e 


= 
= 


1 : 1 


\Ar\ 2 A e , \A e \ 2 A 9 , 

A 2 A—g, A— r A 2 g 


2iD 2 A r — UJ r \K r \A r \ 2 4- Kg \Ag 
(K^ArfAg + K 2 \Ag\ 2 Ag + K 3 A%A_g 


2 ] A r + uj r x 
+ K 4 A- r A 2 g ) 


= 




\Ar\ 2 A r , \M 2 Ag, 
A- r Ag, A 2 A—g 


2iD 2 Ag - UJg [X r \A r \ 2 + \g\Ag 

(LiKI 2 !,. + L 2 \A e \ 2 Ag + L 3 A_ r A 2 


2 ] l e +QJ 9 X 


= 


3 : 1 


A 2 _ x Ag 

- 1 r 


2iD 2 A r - UJ r [K r \A r \ 2 + Kg\Ag\ 2 ] A r - 
2lD 2 Ag-UJg [\ r \A r \ 2 + \g\Ag\ 2 } 


A + ujgLAl 


= 
= 



Table 2: Possible resonances in the third order of approximation. 



oscillating with the frequencies 2uj ri 2uig and \ui r ± uig\. They are given by a particular solution of 
equation 134|) after elimination of secular term and can be expressed as a linear combination 

x 2 = Q^A^A^A^A^g. (36) 

|a|=2 

(2 x) 

Under the assumption of invariancc under the time reflection, constants Q a ' are real and their 

(2 x) — . 

relation to Ca becomes obvious, by substituting x 2 into equation (|34|) . 

When we proceed to the higher order, the discussion is analogical in many respects. The terms 
proportional to e 3 , which appear on the left-hand side of the governing equations, are given by 



8x + uJlx 



a = (Dg + uj 2 )x 3 + 2iu x D 2 A x - 2iuj x D 2 A- x . (37) 



The terms containing Dix± and Dix 2 vanish in consequence of the solvability condition (|35|) . The 
right-hand side contains cubic terms of the Taylor expansion combined using first-order approximations 
r\ , 9\ and quadratic terms combined using one first-order - r\ or 6\ - and one second-order quantity 
- r 2 or 8 2 . Since the second-order terms must be linear combinations of A± r and A±g, the governing 
equations take the form 

{Dl + uj 2 x )x 3 = -2iw x D 2 A x + 2iw x D 2 A_ x + lo 2 x VJ C a 3 ^ A^ Af A a J r A%, (38) 

l«l=3 

where all constants Ca are real. 

The secular terms together with possible resonances are summarized in the table [3 Far from any 
resonance, we eliminate the terms which are secular independently of ui r and uig. Multiplying by 
e -tui x t^ £ ne solvability conditions take the form 

D 2 A r = [n r \A r \ 2 + Kg\A e \ 2 ] A r , (39) 

D 2 Ag = -^[\ r \AA 2 + Xg\Ag\ 2 ]Ag, (40) 
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uj0 : u> r Secular terms Solvability condition 



Outside 
resonance 


D 3 A r 
D 3 Ag 


D 3 A r 

D 3 Ag 


= 
= 


1 : 4 


A 3 

A r Ag 


2iD 3 A r - u) r KA\ 

2iD 3 Ag - LUgLA r Al 


= 
= 


2 : 3 


A r A_ g 

A 2 A 2 


2iD 3 A r - u r KA r A 3 _g 
2iD 3 A e - u e LA 2 r A 2 _g 


= 
= 


3 : 2 


A 2 A 2 
AfA-9 


2iD 3 A r - u r KA 2 _ r A 2 
2iD 3 A g - u e LAlA_g 


= 
= 


4 : 1 


AfAg 


2iD 3 A r - oj r KAlAg 

2iD 3 Ag - OJgLAf 


= 
= 



Table 3: Possible resonances in the fourth order of approximation. 



where we denoted n r = C^q* , Kg — Cnoi > — ^mo an< ^ = C0201 because of simpler notation. 
A particular solution of equation l|38l) is given by linear combination of cubic terms constructed from 
A± r and A±g 

x 3 = E Q { « x) K^^A a J r A a Jg, (41) 

|a|=3 

where all coefficients Qi 3 '^ are real. 

The terms proportional to e 4 in the expanded left-hand side of the equations (|27|) and l|28|) are 



5x + tJix 



= (D 2 + lu 2 x )x 3 + 2D 3 D oXl + 2D Q D 2 x 2 . (42) 



The operator -D0-D2 acts on x-i given by l|3t)l) . The result is found using the solvability conditions l|39|) . 
(14011 and can be written in the form 

2D D 2 x 2 =uj 2 x J2 J^ ] A^AfA a J r A a _%. (43) 

\a\=4 

(x) 

where constants Ja are real because both Do and D 2 produce one i . The right-hand side is 
expanded similarly. We obtain 

(D 2 + w 2 x )x 4 = -2iuj x D 3 A x + 2iu x D 3 A^ x + w 2 x £ C^A^A^A%A a _%, (44) 

l«l=4 

with real C^' x \ On the right-hand side there is only one secular term —2iuj x D 3 A x independently of 
uj r and u>g, the sum contains only terms which becomes secular near a resonance. These terms and 
solvability conditions are listed in the table 

One general feature of an internal resonance k : I is that kw r and liog need not to be infinitesimally 
close. Consider, for example, resonance 1 : 2. The resonance occurs when ujg « 2u r . Suppose that 
the system departs from this exact ratio by small (first-order) deviation tog = 2u r + ea, where a is 
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called detuning parameter. Then the terms A_ r Ag and A 2 in the equations l|34|) remain still secular 
in To since 

A_ r A e = A* r A e i( - Ue - u ^ T ° = A* r A e e i ^ +ta '> T " = A^Age^e^ (45) 
and analogically for A 2 . 

4 The 3:2 parametric resonance 

Let us study oscillations of the conservative system which eigenfrequencies u r and tog are close to 
3 : 2 ratio. The time behavior of the observed frequencies u>* and loq and amplitudes a r and ag of the 
oscillations is given by the solvability conditions l|35|l . <|39|1 and (|40|l . In the fourth order we eliminate 
terms which become nearly secular. For this purpose let us introduce detuning parameters a 2 and 173 
according to 

Zuj r = 2ujg + e 2 a 2 + e 3 a 3 , (46) 

where the term ea\ is missing, because the complex amplitude A depends only on time-scales T 2 and 
T 3 . The secular terms in eq. (I4H) are eliminated if (see table |2J 

2iD 3 A r -u; r a{A 2 r )*A 2 d e- l{ ' T2T2+ ' 73T3) = 0, (47) 
2iD 3 A - cug/3A 3 r A*ge z( ' J2T2+ ' T3T3) = 0, (48) 

where a and [3 are real constants depending on properties of the system. Since A r and Ag are complex, 
the conditions l|47|l and (|48|l together with (|35|) and l|39l) represents 8 real equations. This can be seen 
by substituting the polar forms A r — ^d r e l ^ T and Ag — ^age 1 ^ 6 , and by separating real and imaginary 
parts. We obtain 

D 2 d r = 0, (49) 
D 2 ~ag = 0, (50) 

D 2 (j> r = [n r a 2 . + Kgdg] , (51) 

o 

D 2 <f>8 = ~ [X r a 2 r + Xgd 2 ] , (52) 
o 

D 3 d r = ^^agsin(-3^. + 20fl - <J 2 T 2 - a 3 T 3 ), (53) 

D 3 dg = ^-dldgsm(3<j) r ~2(t)g+a 2 T 2 + c7 3 T 3 ), (54) 
16 

D 3 <f> r = — -~ra r a 2 e cos(-30 r + 2<f>g - <r 2 T 2 - er 3 T 3 ), (55) 
lb 

D 3 (j)g = -^a 3 r cos(3<l> r -2<p + a 2 T 2 +a 3 T 3 ). (56) 

The amplitudes d r and ag of the oscillations change slowly, because they depend only on T 3 . Phases 
4> r and <f>g of oscillations are modified on both time scales T 2 and T 3 . The number of equations can 
be reduced by introducing the phase function 7(12, T 3 ) = 2<j)g — 3(f) r — oT 2 — a 3 T 3 . Then we get 

D 3 a r = — — a r a g smj, (57) 
16 

P^g o 

D 3 dg = —— — a r dg sin7, (58) 
16 

D 2 j = -<r 2 + (fi r aj + \igd 2 g) , (59) 
-D37 = -°"3 + -4-OLr (adg - (3a 2 ) cos 7, (60) 
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were we used the fact that near the resonance ui r ~ (2/3)u>g and we defined fi r — n r — A r and 
l^e = Kg — Xg. The situation can be further simplified if we come back to unique physical time t. Then 
equations for evolution of 7 are merged using d/dt = e 2 D 2 + e 3 D 3 . We obtain 



a g = - 
7 = ' 

where we defined a — ed and a 



—a r a g sin 7, 

P^g 3 
lb 



T 



pi r a 2 + [igc?g + -y (aa 2 e - 13a 2 ) cos 7 



(61) 
(62) 
(63) 



( 2 a 2 + e 3 er 3 . 



4.1 Steady-state solutions 

Steady-state solutions are characterized by constant amplitudes and frequencies of oscillations. Such 
solutions represent singular points of the system governed by equations l|61 [l — 1|63 |) . 

It is obvious from equations l|61|l and (|62|l that the condition d r = dg = can be satisfied (with 
nonzero amplitudes) only if sin 7 = (identically at all times), and thus also 7 = 0. In that case 
equation H63(l transforms to the algebraic equation 



a 

(Jg 



fi r a 2 . + \iga 2 g ± y (aa 2 e - (3a 2 r ) 



The left-hand side can be expressed using the eigenfrequency ratio R — tog/uj r as 



(64) 



R 



(65) 



Then we get 



R= \~^ (Mr«r + ^ea 2 e ) ± ^a r (aa 2 g - [3a 2 ) , 



(66) 



were we neglected terms of order a 4 . Note that the lowest correction to eigenfrequencies is of order 
of a 2 - for given amplitudes a r , ag steady-state oscillations occur when the ratio of eigenfrequencies 
departs from 3/2 by deviation of order of a 2 . 

The relation between observed frequencies of oscillations to*, ojg and eigenfrequencies u> r , ug are 
given by the time derivative of phases <ft r and <pg 



UJ* = LU r + (j) r , UJg — LUg + <j)g. 

We can find simple relation between observed frequencies and the phase function 
3oj* - 2uj* g = 3uj r - 2uj 6 + (34> r - 20 e ) = a + (3<j) r - 2<j) ) = -7. 



(67) 



(68) 



For steady state solutions 7 = 0, and thus observed frequencies are adjusted to exact 3 : 2 ratio even 
if eigenfrequencies depart from it. 

Finally, let us derive explicit relations for oj* and uj up to to the second order in amplitudes. 
Using the equations (|51(l and Ij52(l we find 



1 - - (n r al + KgOg) 



UJg = LJg 



(\ r al + Xga 2 g) 



(69) 
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Figure 2: Co mparison between an analy tical constraint l|71f) and the corresponding numerical solution 
of the system lAbramowicz et all l|2003fl . Each point corresponds to the amplitudes of the oscillations 
at a particular time. On the other hand, from the discussion of equation (|71|) we know that these 
points must lay on an ellipse, whose shape is determined by the multiple-scales method. 



4.2 Integrals of motion 

Behavior of the system is described by three variables a r (t), ag(t) and ^(t) and three first-order 
differential equations (|61|) . I|t)2[) and Ho3f> . However, the number of differential equations can be 
reduced to one because it is possible two find two integrals of motion of the system. Our discussion 
will be analogical to c ase of 1 : 2 resonance of systems with quadratic nonlinearity, as examined by 
iNavfeh & Mookl lll979h . 

Consider equations i|61|) and (|62|l . Eliminating sin 7 from both equations we find 



dt 



(a 2 r + vafj = 



and thus 



va e — const = E, 



(70) 



(71) 



(72) 



where we defined 

au r 2a 

When v > 0, the both amplitudes of oscillations are bounded. The curve [a r (t), ag(t)] is a segment 
of an ellipse. The constant E is proportional to the energy of the system. On the other hand, when 
v < 0, one amplitude of oscillations can grow without bounds while the seco nd amplitude vanishes . 
This case corresponds to the presence of an regenerative element in the system llNavfeh fc Mookll97^ . 
The corresponding curve in the (a r ,ag) plane is a hyperbola. In further discussion we assume that 
v > 0. 

In order to verify that the the energy of the system is conserved, we numerically integrated goy - 
erning equation (|27[) and Q28JI for the one particular system discussed bv lAbramowicz et alJ |2f 
The comparison is in figure [3 The numerical and analytical results are in very good agreement. 
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The second integral of motion is found in following way. Let us multiply the equation i|63|) by ag. 
Then we obtain 

aei = -~oa e + —fi r a r a e + —^ga g + —aa r a e cos 7 —[3a r ag cos 7. (73) 

4 4 o o 

Changing the independent variable from t to ag and multiplying the whole equation by dag we find 

ala 2 9 d(cos-f) + -^-d(a 2 g ) - ^-a 2 r agd{a 2 g ) - ^rd{af) - 
pug p p 

2a 

— —a r a,g cos 'ydag + 2a^ag cos 7 dag = 0. (74) 

P 

The equation (|71|l implies 

agdag = — -. (75) 

v 

With aid of this relation the equation (|74|l takes the form 

ia 2 r a 2 g cos jda r + 2a^a cos jdag + alagd(cos 7) + 

+ ]tg d{al) + % d{at) ~ f d{4) = °" (?6) 

The first three terms express the total differential of function —a^a^ cos 7. Hence, the above equation 
can be arranged to the form 

/ 3 2 8(7 2 , Mr 4 Me 4 

a 1 C ' " 

In other words, 



d I a r a e COS7+ — a e + — a r - — a ) = 0. (77) 



■■) 80" Q LXf A LlQ At — 

a'iag cos 7 + — — a g + —a — a e = const = L (78) 

fjujg 0V P 



is another integral of equations l|61(l - (|63|l . 
4.3 Analytical results 



Knowing two integrals of motion, we should be able to find one differential equation which governs 
behavior of the system. 

First, the amplitudes a r and ag are not independent because they are related by equation (|71|) . 
To satisfy this relation, let us define new variable £(t) by 

a 2 r =£E, a 2 e = (l-^. (79) 

The equation describing an evolution of is derived as follows. Let us multiply equation i|tjlfl by 
2a r and integrate it. We obtain 

d(af) _ a__ 3 2 



dt 

Then we express a 2 using £, and square it. We find 



^) e 2 = («^sin 7 )^. (81) 



uj r a'^ag sin 7. (80) 



au> 7 



The right-hand side of this equation can be expressed using (|78|l 



(a^sin 7 ) 2 = (ala 2 ) 2 - L - ^-a 2 - ^ + ^ajj . (82) 
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Figure 3: The functions ±F(£) = ±(1 — £)£ 3//2 and the quadratic function G(£) which second powers 
are first and second terms on the right-hand side of the equation (|83[) . The behavior of the system 
corresponds to £ in the interval [£i,£a] (denoted by the two dotted vertical lines) where the condition 
> is satisfied. 



After the substitution into the equation 1|81[1 and using the relations l|79|l . we get 

^2 ,, Z?2 



i / m 2 



, e 2 - (i - o 2 e 3 - — 



pvios pv pv z 



(83) 



The equation of motion has very familiar form 

!C 2 £ 2 =F 2 (0-G 2 (0, (84) 

where the K? is a positive constant, F(£) = (1 — £)£ 3//2 an d G(£) is a quadratic function which 
coefficients depend on initial condition through E and L. For example, the equation with an effective 
potential, which governs motion of test particle around a massive body, has the same form. Therefore 
the following discussion is identical as in that case. 

In general, the motion occurs only when £ 2 is positive and thus for £ which satisfy -F 1 ^) > |G(£)|. 
The turning points, where £ changes its signature, are determined by the condition 

\F(>i)\ = |C(OI- (85) 

The functions ±F(£) and G(£) are plotted in figure^ Generally, the function G intersects functions 
±F in two points which corresponds to £(£) oscillating between two bounds £i and £2 given by condition 
l)85[l. In that case the radial and vertical mode of oscillations will periodically exchange the energy. 
The exchanged energy is given by AE/E = £2 — £i- However, for some particular values of L and 
E only one intersection of ±F and G can be found. These stationary oscillations correspond to the 
steady-state solutions discussed above. 

The period of the energy exchange can be find by integration of the equation l|83|l 

T = g (86) 

The integral on the right-hand side can be estimated in the following way. Since P5 = F 2 (£) — G 2 (£) 
is a polynomial of the fifth order in £ having two roots £1 and £2 in the interval [0, 1], we can write it 
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Figure 4: Time evolution of the amplitudes (top panel) and the observed frequencies, v* e — LOg/(2ir) 
(middle panel) and u* — uj*/(2it) (bottom panel). All quantities are rescaled with respect the higher 
eigenfrequency vg . Amplitudes of the oscillations are anticorrelated because the energy is conserved. 
Observed frequencies are correlated because the system is in parametric resonance. 



as — (£ — £i)(£ — £,2)Pz{Qi where P${£) is a polynomial of the third order positive in the interval [0, 1]. 
Using the mean-value theorem we get 

6 di = - f' 2 j = - (87) 

where p > is a value of P3 for some £ in the interval [£1, £2]- Since P5 ~ F 2 ~ 0.01 and (£2 — £i) 2 ~ 
0.01 typically, the values of Pa(£,) are of order of unity, and therefore p ~ 1. The period of the energy 
exchange can be roughly approximated by 

r~^£T 3 / 2 . (88) 

However, near the steady state (£2 ~ £i) 2 approaches to zero and the period becomes much longer. 

The observed frequencies to* and Ug, given by relations (|69() . depend on squares of amplitudes a r 
and ae. Since both a 2 , and a 2 , depend linearly on £(t), also observed frequencies are linear functions 
of £ and are linearly correlated. The slope of this correlation uig = Kui* + Q is independent of the 
energy of oscillations and is given only by parameters of the system, 

iO r K r V — Kg 

The slope of the correlation differs from 3:2, however the observed frequencies are still close to it. 



15 



4.4 Numerical results 



The equations (|fc> 1 1) — (|fc>3|) were solved numerically using the fifth-order Runge-Kutta method with an 
adaptive step size. One of the solutions is shown in figure 0] The top panel of the figure 01 shows the 
time behavior of the amplitudes of two modes of oscillations. Since energy of the system is constant, 
amplitudes are anticorrelated and the two modes are continuously exchanging energy between each 
other. The middle and the bottom panels show the two observed frequencies that are mutually 
correlated. They are also correlated to one of the amplitudes. The frequency ratio varies with time 
and it differs from the exact 3 : 2 ratio, however, it always remains very close to it. The numerical 
solution is in agreement with the general results obtained analytically in the previous section. 



5 Conclusions 

Although this paper was originally motivated by observations and models connected to high-frequency 
QPOs, our results are very general and can be applied to any system with governing equations of the 
form (|27[l and 128[) . Moreover, the solvability conditions, which are derived for all resonances up to 
the fourth order and summarized in tables ^ El an d El are valid also for non-conservative systems. In 
the latter case, the only difference is that constants ci™' 2 ^ that appear in the multiple scale expansion 
are generally complex. However, the results discussed in section^Jwere derived under the assumption 
that the system is conservative and thus all the constants ci" ,a ^ are real. 

The main result of this calculation is the prediction of low-frequency modulation of the amplitudes 
and frequencies of oscill ations. The chara cteristic timescale is approximately given by equation 188|) . 
In a separate paper bv lHorak et alJ (|2JX)4{) we pointed to possible connection of this modulation with 
the 'normal branch oscillations' (NBOs) that are often present together with QPOs. Specifically, we 
suggest that the correlation between the higher frequ ency and t he low er amplitude, evident in figure 
01 is the same as was recently reported in Sco X-l bv lYu et ahl fcnmh. 

It is a pleasure to thank Vladimir Karas, Marek Abramowicz, Wlodek Kluzniak, Paola Rebusco, 
Michal Bursa and Michal Dovciak for helpful discussions. This work was supported by the GACR 
grant 205/03/0902 and GAUK grant 299/2004. 
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